home *** CD-ROM | disk | FTP | other *** search
/ 3D GFX / 3D GFX.iso / amiutils / i_l / irit5 / cagd_lib / bspcoxdb.c < prev    next >
C/C++ Source or Header  |  1995-12-30  |  8KB  |  164 lines
  1. /******************************************************************************
  2. * BspCoxDB.c - Bspline evaluation using Cox - de Boor recursive algorithm.    *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. /*****************************************************************************
  13. * DESCRIPTION:                                                               M
  14. * Returns a pointer to a static data, holding the value of the curve at      M
  15. * the prescribed parametric location t.                                      M
  16. *   Uses the recursive Cox de-Boor algorithm, to evaluate the spline, which  M
  17. * is not very efficient if many evaluations of the same curve are necessary  M
  18. *   Use knot insertion when multiple evaluations are to be performed.        M
  19. *                                                                            *
  20. * PARAMETERS:                                                                M
  21. *   Crv:      To evaluate at the given parametric location t.                M
  22. *   t:        The parameter value at which the curve Crv is to be evaluated. M
  23. *                                                                            *
  24. * RETURN VALUE:                                                              M
  25. *   CagdRType *:  A vector holding all the coefficients of all components    M
  26. *                 of curve Crv's point type. If for example the curve's      M
  27. *                 point type is P2, the W, X, and Y will be saved in the     M
  28. *                 first three locations of the returned vector. The first    M
  29. *                 location (index 0) of the returned vector is reserved for  M
  30. *                 the rational coefficient W and XYZ always starts at second M
  31. *                 location of the returned vector (index 1).                 M
  32. *                                                                            *
  33. * KEYWORDS:                                                                  M
  34. *   BspCrvEvalCoxDeBoor, evaluation, Bsplines                                M
  35. *****************************************************************************/
  36. CagdRType *BspCrvEvalCoxDeBoor(CagdCrvStruct *Crv, CagdRType t)
  37. {
  38.     static CagdRType Pt[CAGD_MAX_PT_COORD];
  39.     CagdBType
  40.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  41.     CagdRType *p, *BasisFunc;
  42.     int i, j, l, IndexFirst,
  43.     k = Crv -> Order,
  44.     Length = Crv -> Length,
  45.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  46.  
  47.     BasisFunc = BspCrvCoxDeBoorBasis(Crv -> KnotVector, k,
  48.                      CAGD_CRV_PT_LST_LEN(Crv),
  49.                      t, &IndexFirst);
  50.  
  51.     /* And finally multiply the basis functions with the control polygon. */
  52.     for (i = IsNotRational; i <= MaxCoord; i++) {
  53.     Pt[i] = 0;
  54.     p = Crv -> Points[i];
  55.     for (j = IndexFirst, l = 0; l < k; )
  56.         Pt[i] += p[j++ % Length] * BasisFunc[l++];
  57.     }
  58.  
  59.     return Pt;
  60. }
  61.  
  62. /******************************************************************************
  63. * DESCRIPTION:                                                               M
  64. * Returns a pointer to a vector of size Order, holding values of the non     M
  65. * zero basis functions of a given curve at given parametric location t.         M
  66. *   This vector SHOULD NOT BE FREED. Although it is dynamically allocated,   M
  67. * the returned pointer does not point to the beginning of this memory and it M
  68. * it be maintained by this routine (i.e. it might be freed next time this    M
  69. * routine is being called).                             M
  70. *   IndexFirst returns the index of first non zero basis function for the    M
  71. * given parameter value t.                             M
  72. *   Uses the recursive Cox de Boor algorithm, to evaluate the Bspline basis  M
  73. * functions.                                     M
  74. *   Algorithm:                                     M
  75. * Use the following recursion relation with B(i,0) == 1.                     M
  76. *                                         M
  77. *          t     -    t(i)            t(i+k)    -   t                        V
  78. * B(i,k) = --------------- B(i,k-1) + --------------- B(i+1,k-1)             V
  79. *          t(i+k-1) - t(i)            t(i+k) - t(i+1)                        V
  80. *                                         M
  81. *   Starting with constant Bspline (k == 0) only one basis function is non   M
  82. * zero and is equal to one. This is the constant Bspline spanning interval   M
  83. * t(i)...t(i+1) such that t(i) <= t < t(i+1). We then raise this constant    M
  84. * Bspline to the prescribed Order and find in this process all the basis     M
  85. * functions that are non zero in t for order Order. Sound simple hah!?         M
  86. *                                                                            *
  87. * PARAMETERS:                                                                M
  88. *   KnotVector:    To evaluate the Bspline Basis functions for this space.   M
  89. *   Order:         Of the geometry.                                          M
  90. *   Len:           Number of control points in the geometry. The length of   M
  91. *                  KnotVector is equal to Len + Order.                       M
  92. *   t:             At which the Bspline basis functions are to be evaluated. M
  93. *   IndexFirst:    Index of the first Bspline basis function that might be   M
  94. *                  non zero.                                                 M
  95. *                                                                            *
  96. * RETURN VALUE:                                                              M
  97. *   CagdRType *:   A vector of length Order thats holds the values of the    M
  98. *                  Bspline basis functions for the given t. A Bspline of     M
  99. *                  order Order might have at most Order non zero basis       M
  100. *                  functions that will hence start at IndexFirst and upto    M
  101. *                  (*IndexFirst + Order - 1).                                M
  102. *                                                                            *
  103. * KEYWORDS:                                                                  M
  104. *   BspCrvCoxDeBoorBasis, evaluation, Bsplines                               M
  105. *****************************************************************************/
  106. CagdRType *BspCrvCoxDeBoorBasis(CagdRType *KnotVector,
  107.                 int Order,
  108.                 int Len,
  109.                 CagdRType t,
  110.                 int *IndexFirst)
  111. {
  112.     static CagdRType
  113.     *B = NULL;
  114.     CagdRType s1, s2, *BasisFunc;
  115.     int i, l, Index,
  116.     KVLen = Order + Len;
  117.  
  118.     if (!BspKnotParamInDomain(KnotVector, Len, Order, FALSE, t))
  119.     CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  120.     if (t == KnotVector[Len])
  121.     t -= IRIT_EPSILON;
  122.     Index = BspKnotLastIndexLE(KnotVector, KVLen, t);
  123.  
  124.     /* Starting the recursion from constant splines - one spline is non     */
  125.     /* zero and is equal to one. This is the spline that starts at Index.   */
  126.     /* As We are going to reference index -1 we increment the buffer by one */
  127.     /* and save 0.0 at index -1. We then initialize the constant spline     */
  128.     /* values - all are zero but the one from t(i) to t(i+1).               */
  129.     if (B != NULL)
  130.     IritFree((VoidPtr) B);
  131.     BasisFunc = B = (CagdRType *) IritMalloc(sizeof(CagdRType) * (1 + Order));
  132.     *BasisFunc++ = 0.0;
  133.     if (Index >= Len + Order - 1) {
  134.     /* We are at the end of the parametric domain and this is open      */
  135.     /* end condition - simply return last point.                */
  136.     for (i = 0; i < Order; i++)
  137.         BasisFunc[i] = (CagdRType) (i == Order - 1);
  138.  
  139.     *IndexFirst = Len - Order;
  140.     return BasisFunc;
  141.     }
  142.     else
  143.     for (i = 0; i < Order; i++)
  144.         BasisFunc[i] = (CagdRType) (i == 0);
  145.  
  146.     /* Here is the tricky part. we raise these constant splines to the      */
  147.     /* required order of the curve Crv for the given parameter t. There are */
  148.     /* at most order non zero function at param. value t. These functions   */
  149.     /* start at Index-order+1 up to Index (order functions overwhole).      */
  150.     for (i = 2; i <= Order; i++) {            /* Goes through all orders... */
  151.     for (l = i - 1; l >= 0; l--) {  /* And all basis funcs. of order i. */
  152.         s1 = (KnotVector[Index + l] - KnotVector[Index + l - i + 1]);
  153.         s1 = APX_EQ(s1, 0.0) ? 0.0 : (t - KnotVector[Index + l - i + 1]) / s1;
  154.         s2 = (KnotVector[Index + l + 1] - KnotVector[Index + l - i + 2]);
  155.         s2 = APX_EQ(s2, 0.0) ? 0.0 : (KnotVector[Index + l + 1] - t) / s2;
  156.  
  157.         BasisFunc[l] = s1 * BasisFunc[l - 1] + s2 * BasisFunc[l];
  158.     }
  159.     }
  160.  
  161.     *IndexFirst = Index - Order + 1;
  162.     return BasisFunc;
  163. }
  164.